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Abstract 



We use a high order accuracy spectral code to carry out two-dimensional 
time-dependent numerical simulations of vortices in accretion disks. In par- 
ticular, we examine the stability and the life time of vortices in circumstellar 
disks around young stellar objects. The results show that cyclonic vortices 
dissipate quickly, while anticyclonic vortices can survive in the flow for hun- 
dreds of orbits. When more than one vortex is present, the anticyclonic 
vortices interact through vorticity waves and merge together to form larger 
vortices. The exponential decay time r of anticyclonic vortices is of the or- 
der of 30 — 60 orbital periods for a viscosity parameter a ~ 10~ 4 (and it 
increases to r m 315 for a = 10~ 5 ), which is sufficiently long to allow heavy 
dust particles to rapidly concentrate in the core of anticyclonic vortices in 
protoplanetary disks. This dust concentration increases the local density of 
centimeter-size grains, thereby favoring the formation of larger scale objects 
which are then capable of efficiently triggering a gravitational instability. 
The relatively long-lived vortices found in this work may therefore play a key 
role in the formation process of giant planets. 

Subject Headings: accretion, accretion disks - circumstellar matter - hy- 
drodynamics - planetary systems - stars: formation - stars: pre-main se- 
quence 
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1 Introduction 



It is presently believed that the formation of the planets in the solar system 
took place via a progressive aggregation of dust grains in the primordial so- 
lar nebula. However, a mechanism for building planetesimals between the 
centimer-size grains (formed by agglomeration and sticking) and the meter- 
size objects capable of triggering a gravitational instability is still lacking 
(see e.g. Adams & Lin 1993 and the recent review of Beckwith, Henning & 
Nakagawa 1999). Recently, it has been suggested (Barge & Sommeria 1995; 
Adams & Watkins 1995; Tanga et al. 1996), that heavy dust particles rapidly 
concentrate in the cores of anticyclonic vortices in the solar nebula, increas- 
ing the local density of centimeter-size grains and favoring the formation of 
larger scale objects which are then capable of efficiently triggering a gravita- 
tional instability. The gravitational instability itself is also enhanced because 
a vortex creates a larger local surface density in the disk. Consequently, even 
if giant planets form (as an alternative model suggests) initially due to a 
gravitational instabilty of the gaseous disk (rather than by planetesimal ag- 
glomeration), vortices may still play an important role. It is the change in 
the Keplerian velocity of the flow in the disk, due to the anticyclonic motion 
in the vortex, that induces a net force towards the center of the vortex (the 
inverse happens in a cyclonic vortex). As a consequence, within a few revo- 
lutions of the vortex around the protostar, the concentration of dust grains 
and the density in the anticyclonic vortex become much larger than outside. 
Analytical estimates (Tanga et al. 1996) have shown that, depending on the 
unknown drag on the dust particles in the gas, the trapping time for dust in 
a small anticyclonic vortex (of size D/r ps 0.01 or less) located at 5 AU (the 
distance of Jupiter from the Sun) is about 20 < r < 10 3 years, i.e. between 
about 2 and 100 local periods of revolution. The mass that can accumulate 
in the core of the vortex on this time scale is of the order of 10~ 5 earth masses 
(i.e. a planetesimal). Barge & Sommeria (1995) considered larger vortices 
(of size D H, where H is the disk half thickness, and H/r ps 0.06 at 5AU) 
and obtained that on a timescale of the order of 500 orbits, a core mass 
of about 16 earth masses can accumulate in the core of the vortex. Since 
small vortices are expected to merge into larger ones, these results should be 
regarded as complementary. The results of Tanga et al. (1996) could repre- 
sent the early stages in planet formation (planetesimals), while the results of 
Barge & Sommeria (1995) could represent a later stage in the evolution (the 
formation of planets cores). 
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It is therefore important for the theory of planet formation to assess un- 
der which conditions vortices form in disks, whether they are stable or not, 
and for how long they can survive in the disk. 

Analytical results (Adams & Watkins 1995, using a geostrophic approxi- 
mation) have shown that many different types of vortex solutions are possible 
in circumstellar disks. However, Adams & Watkins (1995) used simplifying 
assumptions and they were unable to estimate the lifetime and stability of 
the vortices. In order to resolve the issue of the non-linear behavior of vor- 
tices (e.g. vortex merging, scattering, growth, etc.) one needs to carry out 
detailed numerical simulations of the flow. 

recent simplified simulations of a disk (Bracco et al. 1998), which assume 
an incompressible flow and solve the shallow water equations, indicate that 
two-dimensional (large scale) vortices do not fragment into small vortices 
because of the inverse cascade of energy (characteristic of two-dimensional 
flows). Rather, the opposite happens: small vortices merge together to form 
and sustain a large vortex. It is believed for example that the Great Red Spot 
in the atmosphere of Jupiter is a steady solitary vortex also sustained by the 
merging of small vortices (Marcus 1993; see also Ingersoll 1990). However, in 
this case the strong winds in Jupiter's atmosphere are also partially feeding 
vorticity from the background flow into the vortex, and the decay time of 
the vortex is much longer. 

The shear in Keplerian disks inhibits the formation of vortices larger than 
a given scale length L s (L 2 S = v/(dQ/dr), where v is the rotational velocity 
of the vortex (assumed to be subsonic, otherwise the energy is dissipated due 
to sound waves and shocks), and f2 is the angular velocity in the disk). At 
the same time, the viscosity dissipates quickly vortices smaller than a viscous 
scale length L u . Consequently, only vortices of size L satisfying L v < L < L s 
can survive in the flow for many orbits before being dissipated. In the calcula- 
tions of Bracco et al. (1998), positive density perturbation counter rotating 
(anticyclonic) vortices were found to be long lived, while cyclonic vortices 
dissipated very quickly. However, as we noted above, Bracco et al. (1998) 
assumed a shallow water approximation for the disk. These authors do not 
specify the value of the specific heats ratio 7 that was used in the simulations 
(in order for the shallow water equation to represent a polytropic thin disk 
one has to assume 7 = 2; see e.g. Nauta & Toth 1998). Most importantly the 
Mach number of the Keplerian flow (or equivalently the thickness of the disk) 
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is not denned. It is presently not clear whether incompressible results (us- 
ing a shallow water approximation) are valid for more realistic disks, i.e. for 
compressible (and potentially viscous) disks, with a density profile p oc r _15//8 
and a temperature profile T oc r~ 3//4 - the standard Shakura-Sunyaev disk 
model (Shakura & Sunyaev 1973). 

It is the goal of the present work to model vortices in disks, using a 
more realistic approach. We conduct a numerical study of the stability and 
lifetime of vortices in a standard disk model of a protoplanetary nebula. For 
this purpose, we use a time-dependent, two-dimensional high-order accuracy 
hydrodynamical compressible code, assuming a polytropic relation. Here it 
is important to stress the following. The vorticity (circulation) of the flow 
(with a velocity field u) is defined as 

w = Vxu. (1) 

Taking the curl of the Navier-Stokes Equations (see e.g. Tassoul 1978) one 
obtains an equation for the vorticity 

+ u.Vpw oc VP x Vp + ... (2) 

The first term on the right hand side of the equation is the source term for 
the vorticity. Therefore, vorticity can be generated by the non-alignment of 
Vp with VP (e.g. meridional circulation in stars is generated by this term). 
The resultant growth of circulation due to this term is also known as the 
baroclinic instability. However, in the present work we do not address the 
question of the vortex generation process, since we are using a polytropic 
relation P = Kp 1 (K is the polytropic constant, 7 = 1 + 1/n, and n is the 
polytropic index), and the source term vanishes. In this case the vortencity 
(vorticity per unit mass) is conserved (Kelvin's circulation theorem, see e.g. 
Pedlosky 1987) and vortices cannot be generated. 

In the next section we present the numerical modelling of the protoplan- 
etary disk. The results are presented §3 and a discussion follows. 
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2 Accretion Disks modelling 



The time-dependent equations (e.g. Tassoul 1978) were written in cylindri- 
cal coordinates (r,<f>,z), and were further integrated in the vertical (i.e. z) 
direction (e.g. Pringle 1981). We use a Fourier- Chebyshev spectral method 
of collocation (described in Godon 1997) to solve the equations of the disk. 
Spectral Methods are global high-order accuracy methods (Gottlieb & Orszag 
1977; Voigt, Gottlieb & Hussaini 1984; Canuto et al. 1988). These methods 
are fast and accurate and are frequently used to study turbulent flows (e.g. 
She, Jackson & Orszag 1991) and interactions of vortices (Cho & Polvani 
1996a). It is important to stress that spectral codes have very little nu- 
merical dissipation, and that the only dissipation that occurs in the present 
simulations is due to the a viscosity introduced in the equations. All the 
details of the numerical method and the exact form of the equations can be 
found in Godon (1997). We therefore give here only a brief overview of the 
modelling. 

The equations are solved in the plane of the disk (r, 0), with < < 2n 
and R < r < 2R , where the inner radius of the computational domain, R , 
has been normalized to unity. We use an alpha prescription (Shakura & Sun- 
yaev 1973) for the viscosity law, v = ac s H = ac 2 s /Vl K , where c s is the sound 
speed and H = c s /Qk is the vertical height of the disk. The unperturbed 
flow is Keplerian, and we assume a polytropic relation P = Kp 1+1 ^ n (n is the 
polytropic index, while the polytropic constant K is fixed by choosing H/ r at 
the outer radial boundary). Both the inner and the outer radial boundaries 
are treated as free boundaries, i.e. with non-reflective boundary conditions 
where the conditions are imposed on the characteristics of the flow at the 
boundaries. 



The Reynolds number in the flow is given by 

Lu 

tie — , 

V 

where L is a characteristic dimension of the computational domain, u is the 
velocity of the flow (or more precisely, the velocity change in the flow over a 
distance L) and v is the viscosity. Since we are solving for the entire disk, 
L pa r and u ~ vk, and the Reynolds number becomes 

R e = a-\H/r)-\ 
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where we have substituted v = aH 2 Qx, since we are using an a viscosity 
prescription. The Reynolds number in the flow is very high (of the order of 
10 4 — 10 5 for the assumed parameters). 

The simulations were carried out without the use of spectral filters and 
with a moderate resolution of 128 x 128 collocations points. As we noted 
above, the only dissipation in the flow is due to the a viscosity introduced 
in the equations (i.e. the Navier-Stokes equations). 



3 Numerical Results 

In all the models presented here we chose H/r = 0.15 to match protoplan- 
etary disks, however similar results were obtained using H/r = 0.05 and 
H/r = 0.25. The models were first evolved in the radial dimension for an 
initial dynamical relaxation phase (lasting several Keplerian orbits at least). 
Then an axisymmetric disk was constructed from the one-dimensional results, 
on top of which we introduced an initial vorticity perturbation. The initial 
vorticity perturbation had a Gaussian-like form and a constant rotational ve- 
locity of the order of ~ 0.2c s . In all the models we found that cylonic vortices 
dissipate very quickly, within about one local Keplerian orbit, while anticy- 
clonic disturbances persist in the flow for a much longer period of time. It is 
important to stress that in all the models, the anticyclonic vortex, although 
rotating in the retrograde direction (in the rotating frame of reference), has a 
rotation rate that is slower than the local Keplerian flow, and consequently, 
in the inertial frame of reference the vortex rotates in the prograde direction, 
like a planet. As the models are evolved, the initial vorticity perturbation 
is streched and elongated by the shear, within a few Keplerian orbits. The 
elongation of the vorticity into a thin structure transfers enstrophy (the po- 
tential enstrophy is defined as the average of the square of the potential 
vorticity) towards high wave numbers. This process is consistent with the 
direct cascade of enstrophy characteristic of two-dimensional turbulence. It 
forms an elongated negative (relative to the local flow) vorticity strip in the 
direction of the shear, with an elongated vortex in the middle of it (Figure 
1). Due to a Kelvin-Helmoltz instability, perturbations in a forming vortic- 
ity strip propagate (along the strip) in the direction opposite to the shear. 
These propagating waves are called Rossby waves in Geophysics (see e.g. 
Hoskins et al. 1985; in Geophysics this instability is referred to as a shearing 
instability, e.g. Haynes 1987; Marcus, 1993, section 6.2). As a result of this 
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instability the two bendings in the vorticity strip (namely, the trailing and 
the leading vorticity waves of the vortex) move in the direction opposite to 
the shear and fold onto themselves. The extremities of the vorticity waves 
can be again elongated by the shear and and can undergo further folding 
due to the propagation of the Rossby waves. This results in spiral vorticity 
arms around the vortex (Figure 2). The shape of the vortex thus formed, 
then does not change any more during the remaining time of the simulation 
(dozens of orbits). 

We also find that anticyclonic vortices act like overdense regions in the 
disk, i.e. within about one orbit the density increases by about 10 percent in 
the core of the vortex. We cannot check, however, whether a cyclonic vortex 
decreases the density in its core, since cyclonic vortices dissipate very quickly. 



3.1 Flat density profile models 

In a first series of models, we chose a constant density profile throughout 
the disk with a polytropic index n — 3. These models are less realistic 
and are probably closer to the models of Bracco et al. (1998) who used 
an incompressible approximation (the shallow water equation). We ran four 
models with different values of the viscosity parameter a — 1 x 10~ 4 , 3 x 10" 4 , 
6 x 10~ 4 and 1 x 10~ 3 . The viscosity parameter was chosen so as to be 
consistent with values inferred for protostellar disks (e.g. Bell et al. 1995). 
The simulations were followed for up to a maximum time of 60 local Keplerian 
orbits of the vortex in the disk. We found exponential decay times for the 
vortex of r = 3.9 periods for a = 10~ 3 and r = 32.4 periods for a = 10~ 4 
(see Figure 4). In all cases the decay was exponential. 

3.2 Standard disk models 

In an attempt to study the stability and decay times of vortices in more 
realistic disks, we modelled a standard Shakura Sunyaev disk (Shakura & 
Sunyaev, 1973) with p oc r~ 15 / 8 and T cx r~ 3 / 4 (this was achieved by choos- 
ing a polytropic index n = 2.5 together with an ideal gas equation of state). 
In this model we chose a = 10~ 4 , and the simulations were followed for up 
to 20 local Keplerian orbits of the vortex in the disk. Although the initial 
vorticity perturbation was the same as before, the decay time of the vortex 
(r = 60), and its size were found to be larger than for the flat density models 
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(see Figure 4). However, difficulties in assessing precisely the size of the vor- 
tex make it difficult to determine whether the different decay time is due to 
the different size of the vortex or merely to the different density profiles of the 
models. We also ran an additional model with a lower viscosity parameter 
of a — 10~ 5 and found an exponential decay time of more than 315 orbits. 
In this case the resolution was increased to 256 x 256 and the vortex that 
formed was slightly smaller in size than in the previous models. 

In order to gain further insight into the dynamics of vortices, we also ran 
two additional models (one with a = 10~ 4 and one with a = 10~ 3 ) where two 
vorticity perturbations were initially introduced in the flow. In the case of 
a = 10~ 4 , the vortices interact by propagating vorticity waves and eventually 
the two vortices merge together to form a single larger vortex (see Figures 
5a - 5d). In the case of a — 1CT 3 , the vortices do not interact (no vorticity 
waves were observed) and dissipate quickly. 

The results indicate that the exponential decay time of a vortex is in- 
versely proportional to the viscosity (Figure 6), and it can be very large 
indeed (in 30 — 60 orbits the amplitude of the vortex decreases by a factor 
of e) for disks around Young Stellar Objects, where the viscosity might be 
low (a = 10~ 4 ). One expects the decay time to behave like r oc d 2 /u (see also 
Bracco et al. 1998), where d is the size of the vortex and v is the viscosity 
in the flow. In principle, an inviscid model could have vortices that do not 
decay. Therefore, any decay that has previously been observed in inviscid 
numerical simulations (e.g. the shallow water approximation of Bracco et al. 
1998) was probably due to numerical diffusion in the code (the hyperviscosity 
introduced by Bracco et al. 1998 in their model). 



4 Discussion 

Accretion disks possess a very strong shear, which is normally believed to 
lead to a rapid destruction of any structures that form within it (e.g. the 
spiral waves obtained in a perturbed disk by Godon & Livio 1998, dissipate 
or exit the computational domain within about 10 orbits for a = 1 x 10 -4 ). 
However, we find that anticyclonic vortices are surprisingly long-lived, and 
they can survive for hundreds of orbits (the amplitude of the vortex decreases 
exponentially with a time constant of 60 orbits for a = 10~ 4 in disks around 
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young stellar objects). These results are in agreement with other similar sim- 
ulations of the decay of two-dimensional turbulence, e.g. the simulations of 
rotating shallow-water decaying turbulence on the surface of a sphere (Cho 
& Polvani, 1996a & b; modelling of the Jovian atmosphere), where the only 
dissipation of the vorticity is due to the (hyper) viscosity. The size and the 
elongation of the vortices that form out of an initial vorticity perturation in- 
crease with the viscosity. In order for the model to be self consistent, the size 
of the vortices has to be larger than the thickness of the disk (validity of the 
two-dimensional assumption, otherwise the vortices are three dimensional). 
We found that for an alpha viscosity of the order of 1CT 5 (with H/r = 0.15) 
the semi-major axis a of the vortices is slightly smaller than H . However, 
when the viscosity is increased to 10~ 3 the semi-major axis becomes up to 
three times the thickness of the disk. In all the cases the semi-minor axis b 
remains smaller than H, while the elongation (a/b) varies from about 4 (for 
the less viscous case) to about 10 (for the most viscous model). Tanga et 
al. (1996) and Barge & Sommeria (1995) solved numerically the equations 
of motion for dust particles in vortices, and confirmed that dust particles 
concentrate inside vortices on a relatively short timescale. The time taken 
by a dust particle to reach the center of an anticyclonic vortex at a few AU 
ranges from a few orbits to a hundred orbits, depending on the exact value 
of the drag parameter. We have shown that in a standard disk model for a 
protoplanetary disk (a polytropic disk with H/r = 0.15), vortices can survive 
long enough to allow dust particles to concentrate in their core. For a = 10~ 3 
only small vortices would form and would not merge together. In this case 
(using the estimate of Tanga et al. 1996 for small vortices) one finds that the 
vortices would decay within about 10 orbits and the dust concentration in 
the core of the vortices would only reach 10~ 6 Earth masses (planetesimals, 
at 5AU). For a = 10~ 4 small vortices would merge together to form larger 
vortices. The concentration of dust grains in the core of the larger vortices 
could reach about 2 Earth masses (within a hundred revolutions and using 
the estimate of Barge & Sommeria 1995 at 5AU). Therefore, the local density 
of centimeter-size grains could be increased, thus favoring the formation of 
larger scale objects which are then capable of efficiently triggering a gravita- 
tional instability. Our results therefore confirm earlier suspicions that were 
based on a simplified solution of shallow water equations for an incompress- 
ible fluid. 

It is important to note though, that in the present work we have not 
addressed yet the problem of vortex formation. In addition to the baroclinic 
instability described in the Introduction, another potential way to generate 
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vortices in disks around young stellar objects is through infall of rotating 
clumps of gas. It has been suggested that protostellar disks could grow from 
the accretion (or collapse) of rotating gas cloud (e.g. Cassen & Moosman 
1981; Boss & Graham 1993; Graham 1994; Fiebig 1997). The clumps with 
the proper rotation vectors could then give rise to small vortices, that would 
subsequently merge together. 

Finally, we would like to mention that vortices can (in principle at least) 
have many other astrophysical applications. For example, interacting Rossby 
waves can result in radial angular momentum transport (e.g. Llewellyin 
Smith 1996). Vortices are also believed to be important in molecular cloud 
substructure formation in the Galactic disk (e.g. Chantry, Grappin & Leorat 
1993; Sasao 1973). 
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Figures Caption 

Figure 1: A colorscale of the vorticity shows the stretching of an anticy clonic 
vorticity perturbation in a disk model with a viscosity parameter of a = 1CT 4 . This 
results in a negative vorticity strip in the flow, which is then further affected by 
the propagation of Rossby waves along its wings. 

Figure 2: A detailed view of an anticyclonic vortex in the disk, a few orbits 
later than shown in Figure 1 . The vorticity wings have folded up onto themselves 
due to a Kelvin-Helmoltz instability ( accompanied by the propagation of Rossby 
waves). 

Figure 3: The momentum (pv) field is shown in the vicinity of the vortex. The 
coordinates (r,cj)) have been displayed on a Cartesian grid. 

Figure 4-' The maximum amplitude A of the absolute vorticity of the vortex is 
drawn in arbitrary logarithmic units as a function of time (in orbits), for values of 
the alpha viscosity parameter ranging from a = 10~ 5 to a = 10~ 3 . The full lines 
represents the initial flat density profile models, while the doted lines are for an 
initial density profile matching a standard disk model. In all the models, during 
the first orbits, the decay of the maximum amplitude of the vorticity is stronger, 
due to the relaxation of the initial vorticity perturbation. The figure shows that 
the amplitude behaves like A oc e~ l l T , where t, the decay time, increases as the 
viscosity decreases. 

Figures 5a-d. As two vortices interact, they emit vorticity waves and even- 
tually merge to form one large vortex. In this model the viscosity parameter was 
taken to be a = 10 -4 . For convenience the vortices are shown roughly in the same 
orientation in the disk. The complete process of merging takes about 5 — 10 orbits. 

Figure 6: The decay time r against or 1 . The flat-density models of Figure 4 
are represented by stars, together with a straight line r = c/a, where c = 3.2 x 10~ 3 . 
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